close all

%% Data summary
%% theta = 0
% (10,9)
% {'20230827_2000VR_200G_10kHz_K2_Rb2_scan_1'}
% cyclescut = [861 889; 915 937; 1391 1391]
% mol: {{{20230826,17:23},{20230904,23:36},{20230905,2:13}}}
% pkOD: 1065.6 pm 80; molN:11291.3 pm 1160

% (10,10)
% {'20230826_2000VR_200G_10kHz_K2_Rb2_scan_2'}
% cyclescut = []
% mol: {{{20230825,32:41},{20230826,2:16}}}
% pkOD: 1075.4 pm 80; molN: 11308 pm 1248 

%% theta = 0.25 pi
% (10,9)
% {'20230901_2000VR_200G_10kHz_K2_Rb2_scan_1','20230902_2000VR_200G_10kHz_K2_Rb2_scan_2'};
% cyclesCut = [146 151;960 1298]
% molecule: {20230901,40:50},{20230902,2:7},{20230902,34:40}
% pkOD: 1010 pm 68; molN: 10752.8 pm 1190.5

% (10,10)
% {'20230831_2000VR_200G_10kHz_K2_Rb2_scan_1'}
% Cyclescut = [50 63; 71 76; 87 93; 106 111;122 125; 799 853; 1126 1153]
% molecule: {20230831,180:203},{20230901,2:21}
% pkOD: 975.6 pm 53.9; molN: 10599.4 pm 989


%% theta = 0.5 pi
% (10,9)
% {'20230824_2000VR_200G_10kHz_K2_Rb2_scan_2','20230823_2000VR_200G_10kHz_K2_Rb2_scan_1','20230825_2000VR_200G_10kHz_K2_Rb2_scan_1'};
% cyclesCut = [200 1022];
% molecule: {20230823,22:31},{20230824,32:39},{20230825,15:23}
% pkOD: 1113.59 pm 86.6; molN: 11552.4 pm 1251.6 

% (10,10)
% {'20230823_2000VR_200G_10kHz_K2_Rb2_scan_2'}
% Cyclescut = []
% molecule: {20230823,32:40},{20230824,2:16}
% pkOD: 1098.8 pm 65.4; molN: 11697 pm 1049.9


% (9,10)
% {'20230824_2000VR_200G_10kHz_K2_Rb2_scan_1'}
% Cyclescut = []
% molecule: {20230824,21:31}
% pkOD: 1114.9 pm 85.8; molN: 11639.3 pm 1175

% (9,9)
% {'20230822_2000VR_200G_10kHz_K2_Rb2_scan_3'}
% Cyclescut = []
% molecule: {20230822,675:679},{20230823,2:16}
% pkOD: 1063 pm 66; molN: 11100 pm 1133

%% theta = 0.75 pi
% (10,9)
% {'20230902_2000VR_200G_10kHz_K2_Rb2_scan_3','20230903_2000VR_200G_10kHz_K2_Rb2_scan_1'}
% cyclesCut =  [650 1113];
% molecule: {20230902,41:50},{20230903,2:3}
% pkOD: 1057.9 pm 70.8; 11023.9 pm 1210

% (10,10)
% {'20230903_2000VR_200G_10kHz_K2_Rb2_scan_2'}
% cyclesCut = [];
% molecule: {{{20230903,43:51},{20230904,2:16}}}
% pkOD: 1056.5 pm 58.5; molN: 11095.3 pm 959

%% theta = pi
% (10,9)
% {'20230826_2000VR_200G_10kHz_K2_Rb2_scan_1','20230904_2000VR_200G_10kHz_K2_Rb2_scan_2'}
% cyclesCut = [];
% molecule: {{{20230904,23:36},{20230905,2:13},{20230826,17:23}}}
% pkOD: 1065.6 pm 79.8; molN: 11291.3 pm 1160.3

% (10,10)
% {'20230825_2000VR_200G_10kHz_K2_Rb2_scan_2'}
% cyclesCut =  [411 628];
% molecule: {{{20230825,32:41},{20230826,2:16}}}
% pkOD: 1075.4 pm 80.1; molN: 11308 pm 1248



%%
Rot_angle = {0,0,0.25,0.25,0.5,0.5,0.5,0.5,0.75,0.75,1,1}; % [pi]
Files = {{'20230827_2000VR_200G_10kHz_K2_Rb2_scan_1'},{'20230826_2000VR_200G_10kHz_K2_Rb2_scan_2'}, ...
    {'20230901_2000VR_200G_10kHz_K2_Rb2_scan_1','20230902_2000VR_200G_10kHz_K2_Rb2_scan_2'},{'20230831_2000VR_200G_10kHz_K2_Rb2_scan_1'} ...
   {'20230824_2000VR_200G_10kHz_K2_Rb2_scan_2','20230823_2000VR_200G_10kHz_K2_Rb2_scan_1','20230825_2000VR_200G_10kHz_K2_Rb2_scan_1'},  {'20230823_2000VR_200G_10kHz_K2_Rb2_scan_2'},{'20230824_2000VR_200G_10kHz_K2_Rb2_scan_1'}, {'20230822_2000VR_200G_10kHz_K2_Rb2_scan_3'}, ...
   {'20230902_2000VR_200G_10kHz_K2_Rb2_scan_3','20230903_2000VR_200G_10kHz_K2_Rb2_scan_1'}, {'20230903_2000VR_200G_10kHz_K2_Rb2_scan_2'}, ...
   {'20230826_2000VR_200G_10kHz_K2_Rb2_scan_1','20230904_2000VR_200G_10kHz_K2_Rb2_scan_2'}, {'20230825_2000VR_200G_10kHz_K2_Rb2_scan_2'}};
Cycle = {[861 889; 915 937; 1391 1391],[], ...
     [146 151;960 1298],[50 63; 71 76; 87 93; 106 111;122 125; 799 853; 1126 1153], ...
     [200 1022], [], [], [], ...
     [650 1113], [], ...
     [], [411 628] };
Names = {[10,9],[10,10], ...
   [10,9],[10,10], ...
   [10,9],[10,10],[9,10],[9,9], ...
   [10,9],[10,10], ...
   [10,9],[10,10]};

h4 = figure(); 
figure(h4);
% set(h4, 'Position', [scrsz(3) 1 scrsz(3) scrsz(4)])

% in order (10,9), (10,10), (9,10), (9,9)
% coeff = [0.2705 0.3705 0.1984 0.1606]; % considering 532 intensity decay. 
coeff = [0.2709 0.3621 0.2019 0.1651]; % ignore 532 decay

counts = zeros(5,4);

for ss = 1: length(Files)
filename = Files{ss};

%Input parameters
DeltaT = 45*1e3; %units are ns. Represents time between ODT off and REMPI/UV pulse
%DeltaTlist = DeltaT;
DeltaTlist = linspace(0,DeltaT,80);
shotNumStart = 00; % 300
maxShot = 2*10000;

%Changes for decoherence
shotNumEnd = maxShot./2;
Momemtum_filter_scale = 3;
TOF_filter_scale = 3; %3

x_ind=1;

nXmomfilter = 3; %determines how many sigmas to use for filtering criterion. 1 is 1-sigma
nYmomfilter = 3; %determines how many sigmas to use for filtering criterion. 1 is 1-sigma
nTOFfilter = 3; %determines how many sigmas to use for filtering criterion. 1 is 1-sigma
Bfield = 200; %Value of the Bfield in Gauss

state_pair =  Names{ss};
cyclesCut = Cycle{ss};

%%
%Flags for different functions
TOFfilterFlag = 1;
XmomfilterFlag = 1;
YmomfilterFlag = 1;
VariableScanFlag = 0;

CoinShotDistribution = 0;
shot_bin_size = 10000;

%%
TOFAllShots = [];
xPosAllShots = [];
yPosAllShots = [];
shotNumAllShots = [];
xVarAllShots1 = [];
xVarAllShots2 = [];
flagREMPILockAllShots = [];
for i = 1:length(filename)
    filepath = [folderpath, filename{i}];
    if (~exist([filepath,'.dat'],'file'))
        error([filepath,'.dat',' needs to be created using mainMCPtxtAnalysisTiming.m first!'])
    else
        realHit = dlmread([filepath,'.dat']);
    end
    TOFAllShotsSub = realHit(:,1);
    xPosAllShotsSub = realHit(:,2);
    yPosAllShotsSub = realHit(:,3);
    shotNumAllShotsSub = realHit(:,4);
    xVarAllShots1Sub = realHit(:,5);
    xVarAllShots2Sub = realHit(:,6);
    flagREMPILockAllShotsSub = realHit(:,7);
    TOFAllShots = [TOFAllShots;TOFAllShotsSub];
    xPosAllShots = [xPosAllShots;xPosAllShotsSub];
    yPosAllShots = [yPosAllShots;yPosAllShotsSub];
    shotNumAllShots = [shotNumAllShots;shotNumAllShotsSub];
    xVarAllShots1 = [xVarAllShots1;xVarAllShots1Sub];
    xVarAllShots2 = [xVarAllShots2;xVarAllShots2Sub];
    flagREMPILockAllShots = [flagREMPILockAllShots;flagREMPILockAllShotsSub];
end
%%
%Define some relevant constants
mK2amu = 2*(39.96399848); %mass of K2 in amu
mRb2amu = 2*(86.909180520); %mass of Rb2 in amu
A = 15.23; %in mm/sqrt(cm^-1/V)
VR = 2000; %volts


    % Calculate the total number of experimental cycles
    totCycleNum = 1;
    cycleNumAllShots = zeros(length(realHit),1);
    cycleNumAllShots(1) = 1;
    
    for i = 2:length(shotNumAllShots)
        if shotNumAllShots(i) < shotNumAllShots(i-1)
            totCycleNum = totCycleNum + 1;
        end
        cycleNumAllShots(i) = totCycleNum;
    end
    
%     cyclesCut = [1100 1206];[600 1066];[300 2000]; [635 1210]; %[1 429;1160 1190]; % [100 200;300 400]
    [numCutPair,~] = size(cyclesCut);
    % Comments about "cyclesCut"
    % Each pair of values specify an interval of cycles to be cut out of the data analysis
    % Example: to cut from shots 100 - 200 and 300 - 400, use the format cyclesCut = [100 200;300 400];
    % To not cut anything, put [0 0]
    cycleNumInd = 1:length(cycleNumAllShots);


    for i = 1:numCutPair
        tmp = find((cycleNumAllShots <= cyclesCut(i,1))|(cycleNumAllShots >= cyclesCut(i,2))); % find indices of shots that will survive the cyclesCut
        cycleNumInd = intersect(cycleNumInd,tmp);
    end
    % % quick try of bootstrap
    % bootstrap = randi([0,1701],1,800); % cycle number bootstrap index
    % tmp_boot = find(ismember(cycleNumAllShots,bootstrap));
    % cycleNumInd = intersect(cycleNumInd,tmp_boot);
    

%     maxShot = 2*10000;
%     
%     shotNumEnd = maxShot./2; %maxShot./2
    
    %% Define data vectors filtered by shot numbers and cycle numbers
    shotNumIndAll = find((shotNumAllShots >= shotNumStart)&(shotNumAllShots <= shotNumEnd));
    
    filterNumIndAll = intersect(shotNumIndAll, cycleNumInd);
    
    TOFAll = TOFAllShots(filterNumIndAll);
    xPosAll = xPosAllShots(filterNumIndAll);
    yPosAll = yPosAllShots(filterNumIndAll);
    shotNumAll = shotNumAllShots(filterNumIndAll);
    cycleNumAll = cycleNumAllShots(filterNumIndAll);
    xVarAll1 = xVarAllShots1(filterNumIndAll);
    xVarAll2 = xVarAllShots2(filterNumIndAll);
    flagREMPILockAll = flagREMPILockAllShots(filterNumIndAll);
    flagREMPILockAll = ones(length(flagREMPILockAll),1);
    

    % % filter by x-variable 
    VariableValues = unique(xVarAllShots1);
    x_plot = VariableValues(x_ind);
    % filterNumIndAll  = filterNumIndAll(find(xVarAll1 == x_plot));

    LockedIndex = find(flagREMPILockAll == 1 & xVarAll1 == x_plot);
    LockedCycleNumAll = cycleNumAll(LockedIndex);
    LockedCycleNum = length(unique(LockedCycleNumAll));
    
    speciesNames = {'K','K_2','Rb','KRb','K_2Rb','Rb_2','KRb_2','K_2Rb_2'};
    
    speciesMasses = [40,80,87,127,167,174,214,254];

    TOFs  = [11640, 16206, 16870, 20260, 23186, 23656, 26179 - 70, 28400]*1+10400; 
    TOFWindow = 900;
    
    %% ROI window parameters for 99 V and 30 G
    theta = 45;
    if Bfield == 30
        xCentROISpecies = [12.95, 5.024, 5.05, 1.37, 0, -0.72, 0, 0]; %30 Gauss
        yCentROISpecies = [-0.85, -0.5887, -0.85, -0.85, 0, -0.85, 0, 0]; %30 Gauss
    else
        xCentROISpecies = [-11.181, -12.2, -10.9494, -11.333, -11.55, -11.55, -11.55, -11.55]; % 5 Gauss
        yCentROISpecies = [-1.2188, -0.623, -0.5559, -0.8894, -0.8, -0.8, -0.75, -0.75]; % 5 Gauss
    end
    windowROISpecies = [10, 10, 10, 20, 10, 10, 10, 5];
    windowROISpecies = [40, 40, 40, 40, 40, 40, 40, 40]*2;
    xMinROISpecies = xCentROISpecies - windowROISpecies./2;
    xMaxROISpecies = xCentROISpecies + windowROISpecies./2;
    yMinROISpecies = yCentROISpecies - windowROISpecies./2;
    yMaxROISpecies = yCentROISpecies + windowROISpecies./2;
    
    xPosAllRot = cosd(theta)*xPosAll - sind(theta)*yPosAll;
    yPosAllRot = sind(theta)*xPosAll + cosd(theta)*yPosAll;
    
    VariableValues = unique(xVarAllShots1);
    x_plot = VariableValues(x_ind);
    % x_plot

    K2Ind = find((TOFAll >= TOFs(2) - TOFWindow/2) & (TOFAll <= TOFs(2) + TOFWindow/2) ...
        & (xPosAllRot >= xMinROISpecies(2)) & (xPosAllRot <= xMaxROISpecies(2)) ...
        & (yPosAllRot >= yMinROISpecies(2)) & (yPosAllRot <= yMaxROISpecies(2)) ...
        & xVarAll1 == x_plot);
%         & (flagREMPILockAll == 1));
    Rb2Ind = find((TOFAll >= TOFs(6) - TOFWindow/2) & (TOFAll <= TOFs(6) + TOFWindow/2) ...
        & (xPosAllRot >= xMinROISpecies(6)) & (xPosAllRot <= xMaxROISpecies(6)) ...
        & (yPosAllRot >= yMinROISpecies(6)) & (yPosAllRot <= yMaxROISpecies(6)) ...
        & (xVarAll1 == x_plot));
%         & (flagREMPILockAll == 1));
    


    
    TOFK2 = TOFAll(K2Ind);
    xPosK2 = xPosAll(K2Ind);
    yPosK2 = yPosAll(K2Ind);
    shotNumK2 = shotNumAll(K2Ind);
    cycleNumK2 = cycleNumAll(K2Ind);
    meanTOFK2 = mean(TOFK2);
    meanTOFK2Offset = meanTOFK2 + 29944.7;
    
    TOFRb2 = TOFAll(Rb2Ind);
    xPosRb2 = xPosAll(Rb2Ind);
    yPosRb2 = yPosAll(Rb2Ind);
    shotNumRb2 = shotNumAll(Rb2Ind);
    cycleNumRb2 = cycleNumAll(Rb2Ind);
    meanTOFRb2 = mean(TOFRb2);
    meanTOFRb2Offset = meanTOFRb2 + 29944.7;
%     
    COMFitFunction = fittype('a1*(1-exp(-x/a9))+a2','independent','x',...
        'coefficients',{'a1','a2','a9'});
%     COMFitFunction = fittype('a1*x+a2','independent','x',...
%         'coefficients',{'a1','a2'});
        
    
    % exclude the outliers in fitting
    x_mean_K2 = median(xPosK2);
    y_mean_K2 = median(yPosK2);
    x_mean_Rb2 = median(xPosRb2);
    y_mean_Rb2 = median(yPosRb2);
    filter_range =2;
    ind_to_fit_K2_x = find(xPosK2>=x_mean_K2-filter_range & xPosK2<=x_mean_K2+filter_range);
    ind_to_fit_K2_y = find(yPosK2>=y_mean_K2-filter_range & yPosK2<=y_mean_K2+filter_range);
    ind_to_fit_Rb2_x = find(xPosRb2>=x_mean_Rb2-filter_range & xPosRb2<=x_mean_Rb2+filter_range);
    ind_to_fit_Rb2_y = find(yPosRb2>=y_mean_Rb2-filter_range & yPosRb2<=y_mean_Rb2+filter_range);
    
%     COMFitxPosK2 = fit(shotNumK2,xPosK2,COMFitFunction,'Start',[-2.638 5.601 171.6]);
%     COMFitxPosK2 = fit(shotNumK2,xPosK2,COMFitFunction,'Start',[-2.638 5.601]);
    COMFitxPosK2 = fit(shotNumK2(ind_to_fit_K2_x),xPosK2(ind_to_fit_K2_x),COMFitFunction,'Start',[-2.638 5.601 184]);

    
 
    %COMFitxPosK2 = fit(shotNumK2,xPosK2,COMFitFunction,'Start',[-2.638 20 171.6]);
    pointsEvalCOMFitxPosK2 = linspace(shotNumStart,shotNumEnd,shotNumEnd-shotNumStart+1);
    EvalCOMFitxPosK2 = COMFitxPosK2(pointsEvalCOMFitxPosK2);
%     COMFityPosK2 = fit(shotNumK2,yPosK2,COMFitFunction,'Start',[4.688 -8.698 184.7]);
%     COMFityPosK2 = fit(shotNumK2,yPosK2,COMFitFunction,'Start',[4.688 2 184.7]);
%     COMFityPosK2 = fit(shotNumK2,yPosK2,COMFitFunction,'Start',[4.688 2 184.7]);
    COMFityPosK2 = fit(shotNumK2(ind_to_fit_K2_y),yPosK2(ind_to_fit_K2_y),COMFitFunction,'Start',[4.688 2 184.7]);
    pointsEvalCOMFityPosK2 = linspace(shotNumStart,shotNumEnd,shotNumEnd-shotNumStart+1);
    EvalCOMFityPosK2 = COMFityPosK2(pointsEvalCOMFityPosK2);
    RelxPosK2 = xPosK2 - COMFitxPosK2(shotNumK2);
    RelyPosK2 = yPosK2 - COMFityPosK2(shotNumK2);
    
%     COMFitxPosRb2 = fit(shotNumRb2,xPosRb2,COMFitFunction,'Start',[-2.177 1.106 143.9]);
%     COMFitxPosRb2 = fit(shotNumRb2,xPosRb2,COMFitFunction,'Start',[-7.8 0 143.9]);
%     COMFitxPosRb2 = fit(shotNumRb2,xPosRb2,COMFitFunction,'Start',[-7.8 15 143.9]);
    COMFitxPosRb2 = fit(shotNumRb2(ind_to_fit_Rb2_x),xPosRb2(ind_to_fit_Rb2_x),COMFitFunction,'Start',[3 5 184.7]);
    pointsEvalCOMFitxPosRb2 = linspace(shotNumStart,shotNumEnd,shotNumEnd-shotNumStart+1);
    EvalCOMFitxPosRb2 = COMFitxPosRb2(pointsEvalCOMFitxPosRb2);
%     COMFityPosRb2 = fit(shotNumRb2,yPosRb2,COMFitFunction,'Start',[3.407 -3.408 170.7]);
%     COMFityPosRb2 = fit(shotNumRb2,yPosRb2,COMFitFunction,'Start',[3.407 1 170.7]);
    COMFityPosRb2 = fit(shotNumRb2(ind_to_fit_Rb2_y),yPosRb2(ind_to_fit_Rb2_y),COMFitFunction,'Start',[3.407 1 184]);
    pointsEvalCOMFityPosRb2 = linspace(shotNumStart,shotNumEnd,shotNumEnd-shotNumStart+1);
    EvalCOMFityPosRb2 = COMFityPosRb2(pointsEvalCOMFityPosRb2);
    RelxPosRb2 = xPosRb2 - COMFitxPosRb2(shotNumRb2);
    RelyPosRb2 = yPosRb2 - COMFityPosRb2(shotNumRb2);
    
    %%
    shotCycleNumK2 = [cycleNumK2,shotNumK2];
    shotCycleNumRb2 = [cycleNumRb2,shotNumRb2];
    [shotCycleCoin,K2CoinInd,Rb2CoinInd] = intersect(shotCycleNumK2,shotCycleNumRb2,'rows');
    prelimCoinTOFK2 = TOFK2(K2CoinInd);
    prelimCoinshotNumK2 = shotNumK2(K2CoinInd);
    prelimCoincycleNumK2 = cycleNumK2(K2CoinInd);
    prelimCoinRelxPosK2 = RelxPosK2(K2CoinInd);
    prelimCoinRelyPosK2 = RelyPosK2(K2CoinInd);
    prelimCoinTOFRb2 = TOFRb2(Rb2CoinInd);
    prelimCoinshotNumRb2 = shotNumRb2(Rb2CoinInd);
    prelimCoincycleNumRb2 = cycleNumRb2(Rb2CoinInd);
    prelimCoinRelxPosRb2 = RelxPosRb2(Rb2CoinInd);
    prelimCoinRelyPosRb2 = RelyPosRb2(Rb2CoinInd);
    prelimCoinCounts = length(prelimCoinRelxPosK2);
    
    NumShiftPts = 50;
    ShiftVector = -NumShiftPts:1:NumShiftPts;
    EqualShotIndex = ((length(ShiftVector) - 1)/2) + 1;
    ShiftVector = ShiftVector';
    CoinCountsShiftNoFilter = zeros(length(ShiftVector),1);
    CoinCountsShiftMomFilter = zeros(length(ShiftVector),1);
    CoinCountsShiftTOFFilter = zeros(length(ShiftVector),1);
    CoinCountsShiftAllFilter = zeros(length(ShiftVector),1);
    sigmaXmomfilter = 0.41*0.56/Momemtum_filter_scale; %resolution in mm, was 0.41*0.56
    sigmaYmomfilter = 0.41*0.56/Momemtum_filter_scale; %resolution in mm, was 0.41*0.56
    sigmaTOFfilter = 20.5*0.56/TOF_filter_scale; %resolution in ns
    L1 = 0.041;   %[m] length of VMI optics
    L2 = 1 - L1;  %[m] E field free TOF tube length
    eta = ((L2/(2*L1))^2-1);
    RatioK = (eta/meanTOFK2Offset)*DeltaTlist - 1;
    RatioRb = (eta/meanTOFRb2Offset)*DeltaTlist - 1;
    RatioKSquared = RatioK.*RatioK;
    RatioRbSquared = RatioRb.*RatioRb;
    TOFcomparison = nTOFfilter*sigmaTOFfilter*sqrt(RatioKSquared+RatioRbSquared);
    for i = 1:length(ShiftVector)
        shotNumK2Shift = shotNumK2 + ShiftVector(i);
        shotCycleNumK2Shift = [cycleNumK2,shotNumK2Shift];
        [shotCycleCoinShift,K2CoinIndShift,Rb2CoinIndShift] = intersect(shotCycleNumK2Shift,shotCycleNumRb2,'rows');
        prelimCoinTOFK2Shift = TOFK2(K2CoinIndShift);
        prelimCoinRelxPosK2Shift = RelxPosK2(K2CoinIndShift);
        prelimCoinRelyPosK2Shift = RelyPosK2(K2CoinIndShift);
        prelimCoinTOFRb2Shift = TOFRb2(Rb2CoinIndShift);
        prelimCoinRelxPosRb2Shift = RelxPosRb2(Rb2CoinIndShift);
        prelimCoinRelyPosRb2Shift = RelyPosRb2(Rb2CoinIndShift);
        CoinCountsShiftNoFilter(i) = length(prelimCoinRelxPosK2Shift);
        CoinIndXmomfilterShift = find(abs(prelimCoinRelxPosK2Shift + sqrt(mRb2amu/mK2amu)*prelimCoinRelxPosRb2Shift)...
            <= nXmomfilter*sigmaXmomfilter*sqrt(1+(mRb2amu/mK2amu)));
        CoinIndYmomfilterShift = find(abs(prelimCoinRelyPosK2Shift + sqrt(mRb2amu/mK2amu)*prelimCoinRelyPosRb2Shift)...
            <= nYmomfilter*sigmaYmomfilter*sqrt(1+(mRb2amu/mK2amu)));
        MomentumFilterIndexShift = intersect(CoinIndXmomfilterShift,CoinIndYmomfilterShift);
        CoinCountsShiftMomFilter(i) = length(prelimCoinRelxPosK2Shift(MomentumFilterIndexShift));
        prelimCoinRelTOFK2Shift = prelimCoinTOFK2Shift - meanTOFK2;
        prelimCoinRelTOFRb2Shift = prelimCoinTOFRb2Shift - meanTOFRb2;
        WeightedprelimCoinRelTOFK2Shift = prelimCoinRelTOFK2Shift*RatioRb;
        WeightedprelimCoinRelTOFRb2Shift = prelimCoinRelTOFRb2Shift*RatioK;
        TotalprelimCoinRelTOFMatrixShift = abs(WeightedprelimCoinRelTOFRb2Shift + WeightedprelimCoinRelTOFK2Shift);
        TOFDifferenceShift = TotalprelimCoinRelTOFMatrixShift - TOFcomparison;
        TOFConditionalMatrixShift = TOFDifferenceShift <= 0;
        TOFConditionalVectorShift = sum(TOFConditionalMatrixShift,2);
        CoinIndTOFfilterShift = find(TOFConditionalVectorShift > 0);
        CoinCountsShiftTOFFilter(i) = length(prelimCoinRelxPosK2Shift(CoinIndTOFfilterShift));
        FullFilterIndexShift = intersect(CoinIndTOFfilterShift,MomentumFilterIndexShift);
        CoinCountsShiftAllFilter(i) = length(prelimCoinRelxPosK2Shift(FullFilterIndexShift));
        if i == EqualShotIndex
            CoinRelxPosK2 = prelimCoinRelxPosK2Shift(FullFilterIndexShift);
            CoinRelyPosK2 = prelimCoinRelyPosK2Shift(FullFilterIndexShift);
            CoinRelxPosRb2 = prelimCoinRelxPosRb2Shift(FullFilterIndexShift);
            CoinRelyPosRb2 = prelimCoinRelyPosRb2Shift(FullFilterIndexShift);
            prelimCoinCountsindexvector = 1:1:length(prelimCoinRelxPosK2Shift);
            NonCoinCountsIndex = setdiff(prelimCoinCountsindexvector,FullFilterIndexShift);
            NotCoinRelxPosK2 = prelimCoinRelxPosK2Shift(NonCoinCountsIndex);
            NotCoinRelyPosK2 = prelimCoinRelyPosK2Shift(NonCoinCountsIndex);
            NotCoinRelxPosRb2 = prelimCoinRelxPosRb2Shift(NonCoinCountsIndex);
            NotCoinRelyPosRb2 = prelimCoinRelyPosRb2Shift(NonCoinCountsIndex);
        end
    end
    
    BkdCoinCountsShiftNoFilter = CoinCountsShiftNoFilter;
    BkdCoinCountsShiftMomFilter = CoinCountsShiftMomFilter;
    BkdCoinCountsShiftTOFFilter = CoinCountsShiftTOFFilter;
    BkdCoinCountsShiftAllFilter = CoinCountsShiftAllFilter;
    
    BkdCoinCountsShiftNoFilter(EqualShotIndex) = [];
    BkdCoinCountsShiftMomFilter(EqualShotIndex) = [];
    BkdCoinCountsShiftTOFFilter(EqualShotIndex) = [];
    BkdCoinCountsShiftAllFilter(EqualShotIndex) = [];
    
    MeanBkdNoFilter = round(mean(BkdCoinCountsShiftNoFilter),1);
    MeanBkdMomFilter = round(mean(BkdCoinCountsShiftMomFilter),1);
    MeanBkdTOFFilter = round(mean(BkdCoinCountsShiftTOFFilter),1);
    MeanBkdAllFilter = round(mean(BkdCoinCountsShiftAllFilter),1);
    STDBkdNoFilter = round(std(BkdCoinCountsShiftNoFilter),1);
    STDBkdMomFilter = round(std(BkdCoinCountsShiftMomFilter),1);
    STDBkdTOFFilter = round(std(BkdCoinCountsShiftTOFFilter),1);
    STDBkdAllFilter = round(std(BkdCoinCountsShiftAllFilter),1);
    PeakErrorNoFilter = round(sqrt(CoinCountsShiftNoFilter(EqualShotIndex)),1);
    PeakErrorMomFilter = round(sqrt(CoinCountsShiftMomFilter(EqualShotIndex)),1);
    PeakErrorTOFFilter = round(sqrt(CoinCountsShiftTOFFilter(EqualShotIndex)),1);
    PeakErrorAllFilter = round(sqrt(CoinCountsShiftAllFilter(EqualShotIndex)),1);
    DifferenceErrorNoFilter = round(1000*sqrt(CoinCountsShiftNoFilter(EqualShotIndex) + STDBkdNoFilter^2)/LockedCycleNum,1);
    DifferenceErrorMomFilter = round(1000*sqrt(CoinCountsShiftMomFilter(EqualShotIndex) + STDBkdMomFilter^2)/LockedCycleNum,1);
    DifferenceErrorTOFFilter = round(1000*sqrt(CoinCountsShiftTOFFilter(EqualShotIndex) + STDBkdTOFFilter^2)/LockedCycleNum,1);
    DifferenceErrorAllFilter = round(1000*sqrt(CoinCountsShiftAllFilter(EqualShotIndex) + STDBkdAllFilter^2)/LockedCycleNum,1);
   
    

      

    
    %% Plot of coincidence counts after filtering

    subplot(3,4,ss)
    plot(ShiftVector,CoinCountsShiftAllFilter,'b-o','LineWidth',2,'MarkerSize',6)
    hold on
    plot(ShiftVector,MeanBkdAllFilter*ones(length(ShiftVector),1),'r--','LineWidth',1)
    hold off
    xlim([min(ShiftVector), max(ShiftVector)])
    ylim([0, 1.5*max(CoinCountsShiftAllFilter)])
    xlabel('Pulse number difference', 'FontSize', 10)
    ylabel('Ion counts', 'FontSize', 10)
    text(-10, 1.1*max(CoinCountsShiftAllFilter), {['{\it N}_{0} = ',num2str(CoinCountsShiftAllFilter(EqualShotIndex)),' ',char(177),' ',num2str(PeakErrorAllFilter)]}, 'FontSize', 10, 'Color','b')
    %text(0.4*max(ShiftVector), 1.15*max(CoinCountsShiftAllFilter), {['Total Cycles = ',num2str(LockedCycleNum)]}, 'FontSize', 18, 'Color','k')
    text(13, 1.2*max(CoinCountsShiftAllFilter),...
        {['Total Cycles = ',num2str(LockedCycleNum)],'{\it n}_{coin} = {\itn}_{0} - {\it n}_{res}',['          = ',num2str(round(1000*(CoinCountsShiftAllFilter(EqualShotIndex)-MeanBkdAllFilter)/LockedCycleNum,1)),' ',char(177),' ',num2str(DifferenceErrorAllFilter)],...
        '      (cts/1000 cycles)'},...
        'FontSize', 10, 'Color','k')
    text(min(ShiftVector)+2, 0.3*(1.2*max(CoinCountsShiftAllFilter) - MeanBkdAllFilter), {['{\it N}_{res} = ',num2str(MeanBkdAllFilter),' ',char(177),' ',num2str(STDBkdAllFilter)]}, 'FontSize', 10, 'Color','r')
    % text(min(ShiftVector)+2, 1.05*max(CoinCountsShiftAllFilter), {'TOF + Momentum Filtering',['n_{Mom} = ',num2str(nXmomfilter)],['\sigma_{Mom} = ',num2str(sigmaXmomfilter),' mm'],['n_{TOF} = ',num2str(nTOFfilter)],['\sigma_{TOF} = ',num2str(sigmaTOFfilter),' ns']}, 'FontSize', 10, 'Color','k')
    set(gca,'fontsize',10)
    
    set(gcf,'color','white')
    title(['\theta = ',num2str(Rot_angle{ss}),'\pi, State Pair = (',num2str(state_pair(1)),',',num2str(state_pair(2)),')'],'fontweight','bold','FontSize', 10)
    grid off

end



distFig('Screen','Ext')